Interval velocity determination

ABSTRACT

In geophysical exploration, a suite of seismograms is converted to values of interval velocity and dip for each of the subsurface layers. An iteration process fits an arbitrarily dipping Snell&#39;&#39;s Law layered model to the observed field seismograms. Use of the velocity model permits a migration of original data into its true spatial positions.

United States Patent Norman J. Guinzy Austin;

William H. Ruehle, Duncanville, both of Tex.

July 17, 1969 Oct. 5, 1971 t Mobil Oil Corporation Inventors Appl. No.Filed Patented Assignee INTERVAL VELOCITY DETERMINATION 6 Claims, 10Drawing Figs.

US. 340/ 15.5 Int. Cl G0lv 1/00 Field of Search.... 340/155 ReierencesCited UNITED STATES PATENTS 3,417,370 12/1968 Brey PrimaryExaminer-Rodney D. Bennett, Jr.

AssislanLExgminer-H. A. Birmicl Anorneyi-wiiliarn J. Scherback,Frederick E. Dumoulin, Arthur F. Zobalfi Andrew L. Gaboriault and SidneyA. Johnson INPUT T ,AT THRES V N LAYER I TRACE iNOEX E RAY COMPU'PARAMETERS IN FIRST LAYER GENERATE MODEL GENERATE TRAVEL 5 TIME CURVESLOCATE MINJN ERROR CURVE F ERROR THRES.

A N um DROP LAYER R PATENTLEU um 5197;

saw 1 or 5' lNPUT T ,A| ,THRE$.,\ /M N'*LAYER "TRACE INDEX COMPUTE RAY-PARAMETERS IN FIRST LAYER GENERATE MODEL m m N GENERATE TRAVEL TIMECURVES FIG. 1

Loo/kn: MINJN ERROR CURVE ATENTEU [1m 51971 3'611 278 SHEET 2 0F 5 FIG.2

THRES PATENTEDBBT 519m 3511.278

saw u 0F 5 v FIG. 8

FIG. 7

INTERVAL VELOCITY DETERMINATION BACKGROUNDOF THE INVENTION Thisinvention relates to geophysical exploration and more particularly to aprocess for converting field seismograms to the internal velocity andthe dip of each subsurface layer.

In seismic exploration, the measurement of velocity is generallyrecognized as the major parameter in the processing and interpretationof seismograms. The determination of the acoustic velocitycharacteristics from seismograms is described in Seismic Velocities FromSubsurface Measuremerits," C. H. Dix, Geophysics, Vol. 20, pages 68-86,1955.

Briefly, acoustic velocity is usually determined from seismograms by useof the relationship In the foregoing, T,, is the travel time of aseismic wave traveling from a source to an interface and back to asurface receiver. This is the time of the reflection on the seismogram.T is the ideal vertical travel from the surface of the earth to thereflection point. (This is also referred to as'the reflection time atzero horizontal offset.) X is the horizontal offset distance from thesource to the receiver and V, is the acoustic velocity. While the termaverage velocity" is commonly used. it is more accurate to state what isthe RMS acoustic velocity. See Dynamics Correlation Analysis," a paperby W. A. Schneider and Milo Backus, presented at the 36th Annual SEGConvention, Houston, Tex., 1966.

One technique which has recently been used to obtain average velocityfrom a suite of seismograms is described in the Brey US. Pat. No.3,417,370. In a technique such as this, the field data includes a suiteof seismograms each having horizontal offset one with respect to anotheralong a line of exploration. These traces are time shifted one withrespect to another based on different assumed values of averagevelocity. For each reflection, a correlation operation indicates thesignal power between the traces for the different values of assumedaverage velocity. A maximum in the signal power indicates correctaverage velocity to that reflection. This can be carried out foreachreflection on the suite of seismograms so that average velocity toeach of these reflectors is determined.

There are improvements on this type of process in which a plurality ofvalues of dip for each reflector are evaluated. That is, for eachreflector, there is a search through all possible values of dip as wellas velocity to determine the maximum signal power between the traces.This provides outputs representing the dip of each reflector as well asthe average velocity to each reflector. Processes such as this willhereinafter be referred to as a VIP process, an acronym for velocity byintegrated power.

Another processing technique for seismograms includes migration of eachreflection point to its correct depth and horizontal position. Presentlyused migration techniques are employed to account for refraction of theray path in a layered earth. Examples of migration techniques aredescribed in Ex ploration Geophysics, published by Trija Publishing Co.,2nd Edition, I950, Los Angeles, Calif, .I. .I. Jakosky, pages 670 and688. Presently used migration methods employ a velocity which varieswith depth only. Such methods do not properly treat velocity gradientssuch as arise in diverging layers.

SUMMARY OF THE INVENTION The process of the present invention determinesinterval velocities by an iteration process which fits an arbitrarilydipping Snell's Law layered model to the field'data.

Reflections on the seismograms are characterized by three parameters,their zero offset arrival time T their apparent average velocity V andtheir time dip Al. That is, a set of reflections l, 2, 3 i can becharacterized by the set of arrival times T,(X) and dips AT Theseparameters can be determined directly from the field seismograms byseveral processes but the aforementioned VIP-type processes areparticularly suitable for this purpose. The present invention generatesthe interval velocity by fitting a dipping layered model to this data.

In addition, the migrated subsurface attitude of each reflector isobtained. The use of the velocity model permits the migration of theoriginal data into its true spatial position.

DESCRIPTION OF THE DRAWINGS FIG. 1 is a flow sheet depicting the processof the present invention;

FIG. 2 is a field setup suitable for use in obtaining the field data;

FIG. 3a shows ray diagrams; FIG. 3b shows ray diagrams;

FIG. 4 shows the comparison of observed travel times to 1 computedtravel times for one value of velocity;-

' FIG. 5 shows the error function as a function of velocity;

FIG. 6 shows a suite of field seismograms; FIG. 7 depicts data obtainedfrom a VIP-type process; FIG. 8 depicts data obtained from a VIP-typeprocess which includes a dip search; and

FIG. 9 depicts the output of the process of the present invention.

DESCRIPTION OF A PARTICULAR EMBODIMENT exploration. The particular typeof field setup shown in FIG. 2'

is a common-surface-point-recording system. The seismic reflectionsappearing on the suite of seismograms can be approximated by hyperboliccurves across the suite defined by:

X: 7 5"- o n In the foregoing, T is the arrival time of the reflectionon a particular trace in the suite, X is the horizontal offset of theparticular trace, T is the zero offset reflection time and V is averageor apparent velocity.

There are numerous techniques for obtaining apparent average velocityand the reflection time at zero offset. These two parameters are neededto normal moveout correct the suite of seismograms. The VIP-type ofprocess for obtaining apparent velocity and zero offset reflection timewill be briefly described with reference to FIGS. 7 and 8.

A suite of traces similar in appearance to that shown in FIG. 6 isobtained from a field layout. (FIG. 6 includes a number of traces eachhorizontally spaced one from the other. The ordinate in FIG. 6 is recordtime from 0 to about 3.0 seconds. The abscissa represents horizontaloffset along the line of exploration.)

At each increment of record time, the traces are time shifted one withrespect to another. This time shift (or normal moveout correction) isbased upon an assumed value of velocity. Then, the signal power of allof the time-shifted traces is obtained for that value of assumedvelocity. One technique for doing this is to obtain the zero lagcross-correlation function of all traces. The value of the zero lagcross-correlation function is obtained for each assumed value ofvelocity. As shown in FIG. 1.4m: value of average velocity has beeniterated through a plurality of values between 5,000 and 13,000 feet persecond. FIG. 7 shows a plurality of signalpower curves, one signal-powercurve for each of a plurality of record times. Peaks in the signal-powercurve indicate the correct value of assumed apparent velocity for aparticular reflector. In FIG. 7, it will be noted that there are severalreflectors having an apparent velocity indicated to be between 5,000 and7,000 feet per second.

The single curve of FIG. 7 is a plot of the maximum in the signal-pwercurves. This is a good indication of the zero offset time of a pluralityof reflecting layers in the subsurface.

In addition to searching through all values of velocity and verticaltravel time as indicated by the output data of FIG. 7, a process of thistype can also search through a series of assumed values of dip, bothpositive and negative. These components of dip can be expressed eitherin terms of the assumed dipping angle, that is, various values of angle,or they can be expressed in terms of the time shift between the firstand last traces in the suite.

FIG. 8 is a plot of the total energy curves, each similar to the singlecurve of FIG. 7, for each of a plurality of assumed values of clip. Inthis case, the dip has been expressed as total time shift betweenadjacent common depth points. From this information, theapproximate dipof each seismic reflection can be determined.

While the output of a VIP-type process has been depicted pictorially inFIGS. 7 and 8, it will be understood that the outputs will normally bedigitized and available to be used as an input to a process such as thatof the present invention. From a process of the type described above,there is available seismic data including the zero offset arrival timesT of reflections, the apparent velocity V,, for each reflection and thetime dip At, of each seismic reflection. From this data, the process ofthe present invention generates the interval velocity for each layer anda more accurate determination of dip. In addition,

the migrated horizontal and vertical displacement of each reflector isobtained. This is accomplished by generating a. model of the ray pathtaken through the succeeding layers and comparing it to the observedfield data to obtain the best model. This is done layer by layer.

Consider the first layer in conjunction with FIGS. 30 and 3b. Theinterval velocity V of the first layer is determined precisely from:

The average velocity V to the first reflector is known from the outputof the preceding VIP process. This process also specifies the clip 4 (orAr/X) for the first layer. From the dip and the velocity of the firstlayer, the direction of the ray path in the first layer is specified.The zero offset time T from the VIP process specifies the thickness ofthe layer.

Having obtained this information for the first layer directly, now theprocess proceeds to obtain the information for succeeding layers by amodeling procedure. First, different values of velocity for the secondlayer are assumed. That is, the velocity is successively iterated to V VV and so on. For each velocity, Snell's Law specifies the ray directionin the second layer from:

sin a,/ V =sin a,/ V, v

From the ray direction in the second layer, the dip of the second layercan be computed because the ray strikes the second layer normally. Thezero offset time T from the field data, and the iterated intervalvelocity, specifies how thick the second layer Now, the travel time foreach value of X can be computed for each iterated value of velocity.

For example, for the velocity of V the computed reflection times foreach value of horizontal offset are shown by solid line in FIG. 4. Theseare computed from the well-known migration equations; that is, for eachlayer, the horizontal displacement of the ray between the top of theinterval and the bottom of the interval can be determined. The length ofthe ray path in each layer is computed and the travel time over thelayer is determined using the assumed value of interval velocity. 7'

The computed values of travel times T, for different valu of X arerepresented by the solid line in FIG. 4. The observed values ofreflection time on each of the traces are indicated by the circles inFIG. 4. The sum of the differences between the computed and the observedvalues of reflection time for each X are summed to produce an errorfunction. That is, the error Where i is the trace index, m is the numberof traces in the suite, 7', is the observed time of reflection N ontrace 1 and T is the computed time of reflection N on trace i. Thiserror is determined for each of the iterated values of velocity. Therewill be produced a series of values of the error function for each ofthe iterated values of velocity. The error function as a function ofvelocity is shown in FIG. 5. The velocity producing the minimum error isthe proper interval velocity. As will be subsequently explained, theminimum is selected by a Newtonian iteration. The minimum specifies theinterval velocity for the layer.

Having determined the interval velocity V, for the second layer, theprocedure is repeated for the third layer. Summarizing, with referenceto FIG. 3b, the velocities V, and V, are now 'accurately known. The dipof the first and second layers are now accurately known. The zero offsettime T for the third layer indicates the thickness of the third layer.For an assumed value of velocity, the ray direction in the third layeris given by Snell's Law:

Now, the travel times T through the third layer can be computed for eachX. The computation of the travel time through the third layermay be inaccordance with the equation:

Where D; is the length of the ray in 1'" layer. V, is the intervalvelocity in i layer. In computing T the values of interval velocity anddip 0 are used. The travel times for different values of X are computedfor each of a plurality of iterated values of velocity, V V,,, V,,, andso on.

The computed values of T are compared with the observed values of T toproduce an error function. The minimum in this error function selectsthe new value of interval velocity.

At this point, it should be noted that the actual process to besubsequently described is a variation on the foregoing. Specifically,the observed reflection times T from the seismogram are not usedalthough this would be the ideal technique. Rather, their equivalent X 2m a) is used. This term is available for each X as an output of theaforementioned VIP process. It will be seen that the term is related toreflection time by the equation:

However, instead of picking all the reflections from the suite ofseismograms, it is easier to use the average velocity which is availablefrom the VIP process.

The invention may be better understood from the block dia gram ofFIG. 1. The input to the process, indicated at l, includes reflectiontimes T where N is the number of the layer and i is the trace index. Aspointed out previously, instead of reflection times, the input to theprocess may be average velocity from the VIP process. Other inputs toFIG. I include Al-, the dip from the VIP process, THRES, the thresholdvalue for the error function, and V the minimum value of intervalvelocity which will be accepted as being valid.

As indicated at 2, a counter is incremented so that the internalvelocities and dips are successively determined for each layer. Aspreviously indicated, it is quite important that an accuratedetermination of interval velocity for the first layer be obtained. Asindicated at 3, the interval velocity for the first layer is determinedfrom the known average velocity and the dip of the first layer.

As indicated at 4, a model is generated with a successive layer beingadded to the model during each performance of the loop. For each assumedvalue of interval velocity V, Snell's Law specifies the ray direction inthe Nth layer by:

Sin rm-1) This specifies the dip of the Nth layer because the ray pathstrikes the Nth layer normally. 1

Having completed the model down through the Nth layer, the following areknown: The zero offset time T to the Nth layer specifies the thicknessof the layer. The clips of the N-l layer and all previous layers areknown. The interval velocity V,,- of the preceding layer and all of theprevious layers are known. The clip of the Nth layer is known from theSnell's Law computation. Now, the travel times fi] can be computed asindicated at 5. A V

The computed T are subtracted from the corresponding T s from the fielddata to produce an error function. V is iterated through successivevalues to generate an error curve and the minimum in this curve islocated by a Newtonian iteration technique indicated at 6. When theminimum error is located, this error is compared to the thresholdTl-lRES to 'detennine whether it is less than the threshold. This isindicated at 7. lf the error is less than the threshold. then thevelocity mociated with the minimum error, is compared to the minimumacceptab le velocity as indicated at 8. lf 9 is greater thanw then V isstored together with the dip of the Nth layer 0,. as determined from theSnell's Law calculations. The storage step is indicated at 9.

In this case, the counter 2 is incremented and the same functions areperformed for a succeeding layer. If the error threshold comparison at7, or the minimum velocity comparison at 8, indicates that thedetermination of velocity has not been satisfactorily made for thislayer (that is, there is a false indication from the steps 8 or 9), thenthis layer is dropped. This is indicated at 10. For example, if theinput data from the VIP process indicated that there were four layers,1, 2, 3, 4, and there was a false indication in the computation for thethird layer, then this layer is dropped and the layers are numbered sothat the final output will include velocities for layers l, 2 and 3.

From the foregoing, it will be apparent that the method of the presentinvention can be practiced with the use of several well-known types ofcomputing apparatus. The method is particularly suitable for use with ageneral purpose digital computer.

The usefulness of the present invention is best demonstrated withreference to FIGS. 6 through 9.

FIG. 6 is a stacked variable area record section on which the surfacelayers down to 2.0 seconds all dip to the left. At L0 to 1.8 seconds aslight thickening occurs. The reflection labeled Event A" has zero timedip.

Field data from the same region was processed in a VIP- type process toobtain the data depicted in FIGS. 7 and 8. The average velocities andthe T 's from this data were used as input to the process of the presentinvention. The results are displayed in FIG. 9. In FIG. 9, intervalvelocity was determined in accordance with the present invention bothwith dip, as previously described, and with an assumption that there wasno dip. Both of these interval velocities are plotted in FlG.9 as afunction of record time. Also shown in FIG. 9 are arrows which in thisfigure locate the migrated events and give their dip attitude. Note thatEvent A when migrated in spatial coordinates dips to the right thusproducing'a dip reversal not in evidence on the time section of FIG. 6.The data represented by the above figures is summarized in tabular formbelow.

SUMMARY WITH DIP v example, in the foregoing the process was describedin conjunction with layered modeling and Snell's Law ray tracing.Actually, Snell's Law layered'modeling is an approximation, that is, aspecial case, of ray tracing in a continuously varying medium. Ingeneral, Eikonal equations define the ray path in a continuously varyingmedium. These are described, for example, in Mechanical Radiation,"Lindsay, R. B., McGraw Hill, New York 1960, pages 41-47. As applied tothe general case, this invention is applicable to reflecting zonesinstead of layers and the ray paths are determined from the Eikonalequations instead of Snell's Law.

A three-dimensional procws is realized by a simple exten sion of theabove-disclosed process. To the values of T V, and 'b, the value of l,the azimuth, is added. These parameters completely specify the elasticlayering for the threedimensional case, as do the three T V,,, b specifythe twodimensional case. The process is identical in realization to theabove process with the substitution of a three-dimensional modelsubroutine to generate the array of travel times T when N refers to N'"layer, 1 to the f" distance in X direction, j to the f distance in theYdirection, when X and Y are an arbitrary orthogonal coordinate systemon the surface of ground and correspond to the orientation of the arealarray.

While a particular embodiment of the invention has been shown anddescribed, it will be understood that various modifications may be madewithout departing from the true spirit and scope of the invention. Theappended claims are, therefore, intended to cover any suchmodifications.

What is claimed is:

1. The method of accurately determining the interval velocity of eachreflecting zone of a subsurface formation from seismic data representingthe arrival time and time dip of the seismic reflections on tracesobtained along a line of exploration comprising the following steps eachperfonned on automatic computing apparatus:

for a particular reflecting zone N, using the time dip of the Nthseismicreflection to specify the ray path in the upper reflecting zones,determining the ray path in the Nth reflecting zone to determine the dipof the Nth reflecting zone,

generating a computed representation of arrival time for a plurality ofassumed values of interval velocity l and for a plurality of points Xalong said line of exploration,

comparing said computed representation of arrival time with therepresentation of arrival time from said seismic data for each assumedvalue of interval velocity,

picking the value of interval velocity producing the least error in theforegoing step, and

iterating the above steps for successive reflecting zones of i saidformation.

2. The method of accurately determining the interval velocity of eachlayer of a subsurface formation from seismic data representing thearrival time and time dip of the seismic reflections on traces obtainedalong a line of exploration comprising the following steps eachperformed on automatic computing apparatus:

for a particular layer N, using the time dip of the Nth seismicreflection to specify the ray path in the upper layers,

determining from Snell's law at the N-l interface the ray path in theNth layer to determine the dip of 'the Nth interface,

generating a computed representation of arrival time for a plurality ofassumed values of interval velocity V and for a plurality of points Xalong said line of exploration,

comparing said computed representation of arrival time with therepresentation of arrival time from said seismic data for each assumedvalue of interval velocity, picking the value of interval velocityproducing error in the foregoing step, and iterating the above steps forsuccessive layers ofsaid formation.

3. The method recited in claim 2 wherein said representation of arrivaltimcfrom the seismic data include average velocity of each reflectionand wherein average velocity is computed for each assumed value ofinterval velocity and compared with the representations from saidseismic data.

4. The method recited in claim 2 wherein said representations of arrivaltime from'the seismic data include the arrival time of each reflection1),, and wherein reflection time T is computed for each assumed value ofinterval velocity and compared with the representations from saidseismic data.

5. The method recited in claim 2 wherein the step of determining the raypath in the Nth layer is computed in accordance with the following:

the least sin aN=- sin anvelocity determined from a previous iterationof the above steps and 0,, is dip of the N-l interface determined froman iteration of the above steps.

6. The method of accurately determining the interval velocity of eachreflecting zone of a subsurface formation from seismic data representingthe arrival time and time dip of the seismic reflections on tracesobtained along an area of exploration comprising the following stepseach performed on automatic computing apparatus:

for a particular reflecting zone N, using the time dip of the Nthseismic reflection to specify the ray path in the upper reflectingzones, 1

determining the ray path in the Nth reflecting zone to determine the dipand azimuth of the Nth reflecting zone,

generating a computed representation of arrival time for a plurality ofassumed values of interval velocity V,- and for a plurality of points X,Y over said area of exploration.

comparing said computed representation of arrival time with therepresentation of arrival time from said seismic data for each assumedvalue of interval velocity,

picking the value of interval velocity producing the least error in theforegoing step, and

iterating the above steps for successive reflecting zones of saidformation.

poaoso Patent No.

Dated October 5, 1971 Inventor(s) Norman J. Guinzy and William H. RuehleIt is certified that error appears in the above-identified patent andthat said Letters Patent are hereby corrected as shown below:

Column 1, Column 4,

Column 5,

Column 7,

line 6, "internal" should be --interval--.

lines 1-4, the material shown within the quotation marks should bedeleted; and also delete "to"; A

line 6 after "and" change T to T A H II line 27, T should read T lines31-34, that portion of the equation reading T should read T line 32 Vshould read V A line 33 ,"B should read 6 line 25, V should read V H HColumn 8, line 2, B should read 9 A line 16, V should read V Signed andsealed this 28th day of March 1 972.

(SEAL) Attest:

EDWARD M.FLETCHER, JR. Attesting Officer ROBERT GO'ITSCHALK Commissionerof Patents

1. The method of accurately determining the interval velocity of eachreflecting zone of a subsurface formation from seismic data representingthe arrival time and time dip of the seismic reflections on tracesobtained along a line of exploration comprising the following steps eachperformed on automatic computing apparatus: for a particular reflectingzone N, using the time dip of the Nth seismic reflection to specify theray path in the upper reflecting zones, determining the ray path in TheNth reflecting zone to determine the dip of the Nth reflecting zone,generating a computed representation of arrival time for a plurality ofassumed values of interval velocity VN and for a plurality of points Xalong said line of exploration, comparing said computed representationof arrival time with the representation of arrival time from saidseismic data for each assumed value of interval velocity, picking thevalue of interval velocity producing the least error in the foregoingstep, and iterating the above steps for successive reflecting zones ofsaid formation.
 2. The method of accurately determining the intervalvelocity of each layer of a subsurface formation from seismic datarepresenting the arrival time and time dip of the seismic reflections ontraces obtained along a line of exploration comprising the followingsteps each performed on automatic computing apparatus: for a particularlayer N, using the time dip of the Nth seismic reflection to specify theray path in the upper layers, determining from Snell''s law at the N-1interface the ray path in the Nth layer to determine the dip of the Nthinterface, generating a computed representation of arrival time for aplurality of assumed values of interval velocity VN and for a pluralityof points X along said line of exploration, comparing said computedrepresentation of arrival time with the representation of arrival timefrom said seismic data for each assumed value of interval velocity,picking the value of interval velocity producing the least error in theforegoing step, and iterating the above steps for successive layers ofsaid formation.
 3. The method recited in claim 2 wherein saidrepresentation of arrival time from the seismic data include averagevelocity of each reflection and wherein average velocity is computed foreach assumed value of interval velocity and compared with therepresentations from said seismic data.
 4. The method recited in claim 2wherein said representations of arrival time from the seismic datainclude the arrival time of each reflection TNi and wherein reflectiontime TNi is computed for each assumed value of interval velocity andcompared with the representations from said seismic data.
 5. The methodrecited in claim 2 wherein the step of determining the ray path in theNth layer is computed in accordance with the following: where theta N isthe dip of the Nth interface, VN is the assumed value of intervalvelocity for the Nth layer, VN 1 is the interval velocity determinedfrom a previous iteration of the above steps and theta N 1 is dip of theN-1 interface determined from an iteration of the above steps.
 6. Themethod of accurately determining the interval velocity of eachreflecting zone of a subsurface formation from seismic data representingthe arrival time and time dip of the seismic reflections on tracesobtained along an area of exploration comprising the following stepseach performed on automatic computing apparatus: for a particularreflecting zone N, using the time dip of the Nth seismic reflection tospecify the ray path in the upper reflecting zones, determining the raypath in the Nth reflecting zone to determine the dip and azimuth of theNth reflecting zone, generating a computed representation of arrivaltime for a plurality of assumed values of interval velocity VN and for aplurality of points X, Y over said area of exploration, comparing saidcomputed representation of arrival time with the representation ofarrival time from said seismic data for each assumed value of intervalvelocity, picking the value of interval velocity producing the leasterror in thE foregoing step, and iterating the above steps forsuccessive reflecting zones of said formation.